%% test fdtd
% [netfile,path]=uigetfile();
fid=fopen([path,'\',netfile],'r');
nblk=fscanf(fid,'%d',1);
ndata=fscanf(fid,'%d',3);
gg=fscanf(fid,'%f');
gg=reshape(gg,[ndata(:)',3]);
fclose(fid);
x=gg(:,:,:,1);
y=gg(:,:,:,2);
z=gg(:,:,:,3);

nx1=2;nx2=ndata(1)-1;
ny1=2;ny2=ndata(2)-1;
nz1=2;nz2=ndata(3)-1;
delta=10;
% [y,x,z]=meshgrid(delta*(nx1-1:nx2+1),delta*(ny1-1:ny2+1),delta*(nz1-1:nz2+1));
%%
a1x=zeros(size(x)); a1y=zeros(size(x)); a1z=zeros(size(x));
a2x=zeros(size(x)); a2y=zeros(size(x)); a2z=zeros(size(x));
a3x=zeros(size(x)); a3y=zeros(size(x)); a3z=zeros(size(x));
g11=zeros(size(x)); g12=zeros(size(x)); g13=zeros(size(x));
g22=zeros(size(x)); g23=zeros(size(x)); g33=zeros(size(x));
J=zeros(size(x));
%%
a1x(nx1:nx2,ny1:ny2,nz1:nz2)=x(nx1+1:nx2+1,ny1:ny2,nz1:nz2)-x(nx1:nx2,ny1:ny2,nz1:nz2);
a1y(nx1:nx2,ny1:ny2,nz1:nz2)=y(nx1+1:nx2+1,ny1:ny2,nz1:nz2)-y(nx1:nx2,ny1:ny2,nz1:nz2);
a1z(nx1:nx2,ny1:ny2,nz1:nz2)=z(nx1+1:nx2+1,ny1:ny2,nz1:nz2)-z(nx1:nx2,ny1:ny2,nz1:nz2);

a2x(nx1:nx2,ny1:ny2,nz1:nz2)=x(nx1:nx2,ny1+1:ny2+1,nz1:nz2)-x(nx1:nx2,ny1:ny2,nz1:nz2);
a2y(nx1:nx2,ny1:ny2,nz1:nz2)=y(nx1:nx2,ny1+1:ny2+1,nz1:nz2)-y(nx1:nx2,ny1:ny2,nz1:nz2);
a2z(nx1:nx2,ny1:ny2,nz1:nz2)=z(nx1:nx2,ny1+1:ny2+1,nz1:nz2)-z(nx1:nx2,ny1:ny2,nz1:nz2);

a3x(nx1:nx2,ny1:ny2,nz1:nz2)=x(nx1:nx2,ny1:ny2,nz1+1:nz2+1)-x(nx1:nx2,ny1:ny2,nz1:nz2);
a3y(nx1:nx2,ny1:ny2,nz1:nz2)=y(nx1:nx2,ny1:ny2,nz1+1:nz2+1)-y(nx1:nx2,ny1:ny2,nz1:nz2);
a3z(nx1:nx2,ny1:ny2,nz1:nz2)=z(nx1:nx2,ny1:ny2,nz1+1:nz2+1)-z(nx1:nx2,ny1:ny2,nz1:nz2);

g11=a1x.*a1x+a1y.*a1y+a1z.*a1z;
g12=a1x.*a2x+a1y.*a2y+a1z.*a2z;
g13=a1x.*a3x+a1y.*a3y+a1z.*a3z;
g22=a2x.*a2x+a2y.*a2y+a2z.*a2z;
g23=a2x.*a3x+a2y.*a3y+a2z.*a3z;
g33=a3x.*a3x+a3y.*a3y+a3z.*a3z;

J=a1x.*(a2y.*a3z-a2z.*a3y)+a2x.*(a3y.*a1z-a3z.*a1y)+a3x.*(a1y.*a2z-a1z.*a2y);

%%
ex=zeros(nx2-nx1+3,ny2-ny1+3,nz2-nz1+3);ey=ex;ez=ex;
cex=ex;cey=ex;cez=ex;chx=ex;chy=ex;chz=ex;
rand('seed',19860217);
chx=rand(nx2-nx1+3,ny2-ny1+3,nz2-nz1+3)*0;
chy=rand(nx2-nx1+3,ny2-ny1+3,nz2-nz1+3)*0;
chz=rand(nx2-nx1+3,ny2-ny1+3,nz2-nz1+3)*0;
chz(floor(nx1+nx2)/2,floor(ny1+ny2)/2,2)=100;
chx(:,:,1)=0;chx(:,:,nz2+1)=0;chx(1,:,:)=0;chx(nx2+1,:,:)=0;chx(:,1,:)=0;chx(:,ny2+1,:)=0;
chy(:,:,1)=0;chy(:,:,nz2+1)=0;chy(1,:,:)=0;chy(nx2+1,:,:)=0;chy(:,1,:)=0;chy(:,ny2+1,:)=0;
chz(:,:,1)=0;chz(:,:,nz2+1)=0;chz(1,:,:)=0;chz(nx2+1,:,:)=0;chz(:,1,:)=0;chz(:,ny2+1,:)=0;
%%
sigmaex=0.01;sigmaey=sigmaex;sigmaez=sigmaex;
mu=4*pi*1E-7;dt=1E-6;gamma=10.0/mu*(dt/delta)^2;
tmax=1E-3;
%%
pp=[];
for n=1:floor(tmax/dt)
    hx((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g11((nx1:nx2),(ny1:ny2),(nz1:nz2)).*chx((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g12((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chy((nx1:nx2)-1,(ny1:ny2)+1,(nz1:nz2)) ...
        +chy((nx1:nx2),(ny1:ny2)+1,(nz1:nz2))...
        +chy((nx1:nx2)-1,(ny1:ny2),(nz1:nz2))...
        +chy((nx1:nx2),(ny1:ny2),(nz1:nz2))    )...
        +1/4*g13((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chz((nx1:nx2)-1,(ny1:ny2),(nz1:nz2)+1) ...
        +chz((nx1:nx2),(ny1:ny2),(nz1:nz2)+1)...
        +chz((nx1:nx2)-1,(ny1:ny2),(nz1:nz2))...
        +chz((nx1:nx2),(ny1:ny2),(nz1:nz2)));

    hy((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g22((nx1:nx2),(ny1:ny2),(nz1:nz2)).*chy((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g23((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chz((nx1:nx2),(ny1:ny2)-1,(nz1:nz2)) ...
        +chz((nx1:nx2),(ny1:ny2),(nz1:nz2)+1)...
        +chz((nx1:nx2),(ny1:ny2)-1,(nz1:nz2)+1)...
        +chz((nx1:nx2),(ny1:ny2),(nz1:nz2))    )...
        +1/4*g12((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chx((nx1:nx2),(ny1:ny2)-1,(nz1:nz2)) ...
        +chx((nx1:nx2)+1,(ny1:ny2),(nz1:nz2))...
        +chx((nx1:nx2)+1,(ny1:ny2)-1,(nz1:nz2))...
        +chx((nx1:nx2),(ny1:ny2),(nz1:nz2)));

    hz((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g33((nx1:nx2),(ny1:ny2),(nz1:nz2)).*chz((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g23((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chy((nx1:nx2),(ny1:ny2),(nz1:nz2)-1) ...
        +chy((nx1:nx2),(ny1:ny2)+1,(nz1:nz2))...
        +chy((nx1:nx2),(ny1:ny2)+1,(nz1:nz2)-1)...
        +chy((nx1:nx2),(ny1:ny2),(nz1:nz2))    )...
        +1/4*g13((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        chx((nx1:nx2),(ny1:ny2),(nz1:nz2)-1) ...
        +chx((nx1:nx2)+1,(ny1:ny2),(nz1:nz2))...
        +chx((nx1:nx2)+1,(ny1:ny2),(nz1:nz2)-1)...
        +chx((nx1:nx2),(ny1:ny2),(nz1:nz2)));

    cex((nx1:nx2),(ny1:ny2),(nz1:nz2)) = 2*dt*( ...
        + hy( (nx1:nx2),(ny1:ny2),-1 + (nz1:nz2)) ...
        - hy( (nx1:nx2),(ny1:ny2),(nz1:nz2)     ) ...
        - hz( (nx1:nx2),-1 + (ny1:ny2),(nz1:nz2)) ...
        + hz( (nx1:nx2),(ny1:ny2),(nz1:nz2)     ) ...
        )./(J((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma + dt*sigmaex))...
        + (cex((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma - dt*sigmaex))/(2*gamma + dt*sigmaex);

    cey((nx1:nx2),(ny1:ny2),(nz1:nz2)) = 2*dt*(...
        - hx( (nx1:nx2),(ny1:ny2),-1 + (nz1:nz2))...
        + hx( (nx1:nx2),(ny1:ny2),(nz1:nz2)     )...
        + hz( -1 + (nx1:nx2),(ny1:ny2),(nz1:nz2))...
        - hz( (nx1:nx2),(ny1:ny2),(nz1:nz2))...
        )./(J((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma + dt*sigmaey))...
        +  (cey((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma - dt*sigmaey))/(2*gamma + dt*sigmaey);

    cez((nx1:nx2),(ny1:ny2),(nz1:nz2)) = 2*dt*(...
        + hx((nx1:nx2),-1 + (ny1:ny2),(nz1:nz2)) ...
        - hx((nx1:nx2),(ny1:ny2),(nz1:nz2)) ...
        - hy(-1 + (nx1:nx2),(ny1:ny2),(nz1:nz2)) ...
        + hy((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        )./(J((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma + dt*sigmaez))...
        +  (cez((nx1:nx2),(ny1:ny2),(nz1:nz2))*(2*gamma - dt*sigmaez))/(2*gamma + dt*sigmaez);

    ex((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g11((nx1:nx2),(ny1:ny2),(nz1:nz2)).*cex((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g12((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cey((nx1:nx2),(ny1:ny2),(nz1:nz2)) ...
        +cey((nx1:nx2),(ny1:ny2)-1,(nz1:nz2))...
        +cey((nx1:nx2)+1,(ny1:ny2)-1,(nz1:nz2))...
        +cey((nx1:nx2)+1,(ny1:ny2),(nz1:nz2))    )...
        +1/4*g13((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cez((nx1:nx2)+1,(ny1:ny2),(nz1:nz2)) ...
        +cez((nx1:nx2)+1,(ny1:ny2),(nz1:nz2)-1)...
        +cez((nx1:nx2),(ny1:ny2),(nz1:nz2)-1)...
        +cez((nx1:nx2),(ny1:ny2),(nz1:nz2)));

    ey((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g22((nx1:nx2),(ny1:ny2),(nz1:nz2)).*cey((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g23((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cez((nx1:nx2),(ny1:ny2)+1,(nz1:nz2)) ...
        +cez((nx1:nx2),(ny1:ny2),(nz1:nz2)-1)...
        +cez((nx1:nx2),(ny1:ny2)+1,(nz1:nz2)-1)...
        +cez((nx1:nx2),(ny1:ny2),(nz1:nz2))    )...
        +1/4*g12((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cex((nx1:nx2),(ny1:ny2)+1,(nz1:nz2)) ...
        +cex((nx1:nx2)-1,(ny1:ny2),(nz1:nz2))...
        +cex((nx1:nx2)-1,(ny1:ny2)+1,(nz1:nz2))...
        +cex((nx1:nx2),(ny1:ny2),(nz1:nz2)));

    ez((nx1:nx2),(ny1:ny2),(nz1:nz2)) = g33((nx1:nx2),(ny1:ny2),(nz1:nz2)).*cez((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        +1/4*g23((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cey((nx1:nx2),(ny1:ny2),(nz1:nz2)+1) ...
        +cey((nx1:nx2),(ny1:ny2)-1,(nz1:nz2))...
        +cey((nx1:nx2),(ny1:ny2)-1,(nz1:nz2)+1)...
        +cey((nx1:nx2),(ny1:ny2),(nz1:nz2))    )...
        +1/4*g13((nx1:nx2),(ny1:ny2),(nz1:nz2)).* (...
        cex((nx1:nx2),(ny1:ny2),(nz1:nz2)+1) ...
        +cex((nx1:nx2)-1,(ny1:ny2),(nz1:nz2))...
        +cex((nx1:nx2)-1,(ny1:ny2),(nz1:nz2)+1)...
        +cex((nx1:nx2),(ny1:ny2),(nz1:nz2)));
    pp(n)=max(abs(cex(:)));
    chx((nx1:nx2),(ny1:ny2),(nz1:nz2)) = ( ...
        -(dt*ey((nx1:nx2),(ny1:ny2),(nz1:nz2)))/(mu)...
        + (dt*ey((nx1:nx2),(ny1:ny2),1 + (nz1:nz2)))/(mu) ...
        + (dt*ez((nx1:nx2),(ny1:ny2),(nz1:nz2)))/(mu) ...
        - (dt*ez((nx1:nx2),1 + (ny1:ny2),(nz1:nz2)))/(mu) ...
        )./J((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        + chx((nx1:nx2),(ny1:ny2),(nz1:nz2));

    chy((nx1:nx2),(ny1:ny2),(nz1:nz2)) = (...
        (dt*ex((nx1:nx2),(ny1:ny2),(nz1:nz2)))/(mu)...
        - (dt*ex((nx1:nx2),(ny1:ny2),1 + (nz1:nz2)))/(mu)...
        - (dt*ez((nx1:nx2),(ny1:ny2),(nz1:nz2)))/(mu) ...
        + (dt*ez(1 + (nx1:nx2),(ny1:ny2),(nz1:nz2)))/(mu) ...
        )./J((nx1:nx2),(ny1:ny2),(nz1:nz2))...
        + chy((nx1:nx2),(ny1:ny2),(nz1:nz2));
    chz(:,:,nz2+1)=0;
    for k=nz2:-1:nz1
        chz((nx1:nx2),(ny1:ny2),k) = -chx((nx1:nx2),(ny1:ny2),k)...
            + chx(1 + (nx1:nx2),(ny1:ny2),k) ...
            - chy((nx1:nx2),(ny1:ny2),k) ...
            + chy((nx1:nx2),1 + (ny1:ny2),k) ...
            + chz((nx1:nx2),(ny1:ny2),1 + k);
    end

    surf(cex(nx1:nx2,ny1:ny2,2));
    drawnow;


end


